Population genetics and evolutionary history of the intertidal brittle star Ophiothrix (Ophiothrix) exigua in the northern China Sea

Abstract Ophiothrix (Ophiothrix) exigua is a common brittle star in the northwestern Pacific. As a dominant species, O. exigua inhabiting the intertidal rocky ecosystem are affected by multiple environmental stressors, but molecular insights into their genetic population structure remain poorly studied. In this study, we investigated the population genetics and evolutionary history of six O. exigua populations from the northern China Sea using mitochondrial (COI, NAD4) and nuclear (ITS2, 18S) gene markers. High haplotype diversity, low nucleotide diversity, and low rates of gene differentiation among the populations of O. exigua were detected. Pairwise genetic differentiation (ΦST) statistics between different localities were negative or low and insignificant, suggesting strong gene flow of this species over the study areas. The phylogenetic analyses showed that the populations exhibited high homogeneity between localities in our study area. Demographic analyses indicated that the populations experienced sustained expansion around 0.2 million years ago. This expansion was likely related to transgressions events in the Yellow Sea during the Pleistocene period. Additional samples of O. exigua from disparate geographical locations, especially the Japan Sea and the Korean Peninsula, will be needed to unravel the population genetic patterns and evolutionary history of this species.


| INTRODUC TI ON
Genetic diversity is a vital component of biodiversity and a key parameter for assessing the status of biological resources (Cen et al., 2023).Rich genetic diversity helps enhance the adaptability and evolutionary potential of species (Hughes et al., 2008;Yang et al., 2016).Various processes such as divergence, isolation, and bottleneck all influence the levels of genetic diversity of population (Raman et al., 2023).However, the genetic diversity, population structure, and ecological functions of marine organisms have all been profoundly affected by global warming and human activities (Donelson et al., 2019;Eirin-Lopez & Putnam, 2019;Törnroos et al., 2019).Rising ocean temperatures pose a significant challenge to the survival of marine organisms (Schiel et al., 2004), potentially leading to the northward migration of some benthic organisms (Xu et al., 2023).Assessing and maintaining the genetic diversity of marine organisms is crucial for preserving biodiversity and sustaining the health of marine ecosystems.
The Yellow and Bohai Seas, as important continental marginal seas in the northwestern Pacific Ocean (Liu, 2013), are sensitive to global change due to the wide, flat, and low-gradient shelves (Liu et al., 2016(Liu et al., , 2018)).Since the late Quaternary, the frequent "transgression" and "regression" events in the Yellow and Bohai Seas have been primarily driven by the significant fluctuations in sea level caused by climate change (Liu et al., 2014;Yao et al., 2017).
Sedimentary environment studies in the Yellow and Bohai Sea region have revealed that this area has experienced at least four marine transgression events since 200 Ka (Yao et al., 2022), which likely correspond to the interglacial periods indicated by Marine Isotope Stages (MIS)7, 5, 3, and 1 (Yao et al., 2022;Yi et al., 2013).These cyclic transgression and regression events have significantly impacted the population structure of marine species inhabiting the Yellow and Bohai Sea Shelf (Avise, 2000;Han et al., 2008;Shen et al., 2011;Yang & Li, 2016), leading to population growth and range expansion (Han et al., 2008;Wang, 1999) during transgression or population bottlenecks, genetic drift, and fragmented or confined to refugia during regression (Liu et al., 2006;Voris, 2000).The alternation between transgression and regression events notably affected intertidal species, likely shaping their current distribution patterns and driving their rapid adaptation to ever-changing environments (Xue et al., 2014).
In this study, we aim to elucidate the population genetics and evolutionary history of O. exigua from the northern China Sea using DNA sequence data from both mitochondrial and nuclear genes to (1) assess the genetic structure among different geographic populations; (2) examine the genetic diversity; and (3) infer the demographic changes of O. exigua.

| Sampling and DNA extraction
A total of 80 samples of O. exigua were collected from six intertidal zones of locations including Qinhuangdao (QHD), Dalian (DL), Rongcheng (RC), Langyatai (LYT), Renjiatai (RJT), and May Fourth Square (WS, Figure 1).The morphology of O. exigua exhibits slight variation between different sampling points (Table S1).After morphological identification, the collected samples were fixed with absolute ethanol and stored at −20°C.A segment of approximately 1 cm from the arm tissue of the brittle star was taken for DNA extraction.After grinding in liquid nitrogen, total genomic DNA was manuals.The extracted total genomic DNA was qualified by 1% agarose gel electrophoresis before PCR amplification.

| PCR amplification and sequencing
For each sample, we conducted amplification with specific primers for two mitochondrial genes-cytochrome c oxidase subunit I (COI) and NADH dehydrogenase subunits 4 (NAD4)-and two nuclear genes-Internal Transcribed Spacer 2 (ITS2) and 18S ribosomal RNA (18S) (Table 1).The primers for the NAD4 were designed based on the complete mitochondrial genome sequence of O. exigua (GenBank number: ON729296) (Table 1).PCRs were conducted in a total volume of 25 μL mixture, comprising 1 μL of each primer (10 μM), 9.5 μL of diluted distilled H 2 O, 12.5 μL of 2× Taq PCR Master Mix (with blue dye, Sangon Biotech, Inc), and 1 μL of the template DNA.The PCR cycling conditions for COI, ITS2, and 18S were followed as described in the respective references (Table 1).The amplification of NAD4 used the following profile: The first step was at 95°C for 3 min, followed by 35 cycles at 95°C for 30 s, 50°C for 30 s and 72°C for 1 min, and a further extension at 72°C for 10 min.The PCR products were sequenced for both directions by Sangon Biotech, Inc. (Qingdao, China) with the same primers used in the amplification reaction.All sequences of this study were uploaded to the GenBank database (Table 2).
The alignment of each gene fragment was then auto-trimmed by G-blocks (Jose, 2000) with default settings to output conserved datasets.Genetic diversity indices, including the number of nucleotide composition, segregating sites (S), number of haplotypes (H), nucleotide diversity (π), haplotype diversity (Hd), and gene flow, were estimated for each dataset using DnaSP v6.12.03 (Librado & Rozas, 2009).The 18S dataset showed highly conserved intraspecific variation, with populations from the six locations sharing a single haplotype.Due to the lack of informative variation, the 18S marker was excluded from further analyses.
The Φ ST values were estimated using a pairwise difference distance matrix, and significance levels were assessed by 1000 permutations.

| Genetic variation
We calculated genetic diversity statistics for the six locations of O. exigua using three gene fragments (Table 3).The populations presented high haplotype diversity, with Hd values ranging from 0.995 to 0.998, and low nucleotide diversity, with π values between 0.00882 and 0.01372.NAD4 exhibited the highest nucleotide diversity (π = 0.01372), followed by COI (π = 0.00985).In contrast, ITS2 appeared to be the most conserved gene in O. exigua compared with the other two fragments (π = 0.00032).Six populations all presented high genetic diversity according to their high haplotype diversities (Hd: 0.933-1).The nucleotide diversities ranged from 0.00760 to 0.01622, with the LYT population exhibiting the highest values for both single mitochondrial gene (π = 0.01119 for COI, π = 0.01622 for NAD4) and the concatenated genes (π = 0.01074).

| Phylogenetic and network analyses
The median-joining network (Figure 2 To study the phylogenetic relationship of O. exigua in six regions, we constructed the haplotype phylogenetic trees using the concatenated COI-NAD4-ITS2 gene set for all populations.ML and BI analyses produced phylogenetic trees with consistent topology but with different supporting values (Figure S1).The phylogenetic analysis of O. exigua populations revealed a lack of stable geographical differentiation.
Despite the presence of unique haplotypes within each population, individuals from different locations were intermixed in the phylogenetic tree, without any obvious group clustering observed.

| Genetic structure
The AMOVA results showed that a large proportion of molecular genetic variation was present within populations (99.25%, p < .001),and the remaining variation was found among populations (0.75%, p < .001)for the concatenated sequences.The overall fixation index (F ST ) observed for all populations was nonsignificant (0.00753, p > .05,Table 4).Consistent with the concatenated dataset, the genetic variance percentages were high for all three markers: 99.94% for COI, 100.12% for NAD4, and 99.25% for ITS2 (Table 5).Pairwise Φ ST values between populations varied from −0.009 (QHD-WS) to

TA B L E 2
The information of GenBank accession.

Locus GenBank accession number
0.081 (DL-RJT) for the concatenated dataset.Population pairwise Φ ST analysis among all populations demonstrated no significant genetic differentiation.The mean K2P distance was between 0.008 and 0.010 for all populations (Table 5).

| Historical dynamics of the population
BSP analysis, mismatch distribution, and neutrality tests were employed to investigate the historical demography of the O. exigua population.The overall negative neutrality tests (Tajima's D and Fu's Fs) suggest that O. exigua encountered a population bottleneck followed by expansion (Table 3).The mismatch distribution curves drawn based on COI and NAD4 mitochondrial gene sequences showed unimodal distribution, which was consistent with the population expansion model (Figure 3).Using the formula t = 2μk/τ (τ = 5.6), the time of the expansion event was calculated as 193,327 years before present (BP).BSP analysis of the total population also showed an increase in the effective population size, indicating that the populations experienced sustained expansion at about 0.2 Ma (Figure 4) and that the population rapidly grew to its current size.

| DISCUSS ION
In the present study, we found high population homogeneity across the Yellow and Bohai Seas.The populations of O. exigua in the studied areas exhibit a relatively high level of genetic diversity, but no population divergence (Figure 2 and Table 3).Our results demonstrated that the different morphs present in O. exigua correspond to a single evolutionary unit, with no evidence of cryptic species, as we observed no genetic divergence was observed between those variant individuals.The homogeneity in O. exigua contrast with the high number of cryptic speciation events detected in other ophiothrix, such as Ophiothrix fragilis f. nuda Madsen, 1970(Pérez-Portela et al., 2013;Taboada & Pérez-Portela, 2016) and Ophiothrix angulata Say, 1825(Hernández-Díaz et al., 2023).The difference in morphological variation of O. exigua may be the result of phenotypic plasticity due to adaptation to the environment.Furthermore, the extremely high genetic diversity found in O. exigua seems to be a characteristic of some abundant ophiuroids likely related to large population sizes (Hunter & Halanych, 2010), which makes populations less susceptible to genetic drift and helps maintain a high-level of low-frequency haplotypes within populations (Muths et al., 2009;Pérez-Portela et al., 2013).A higher genetic diversity also helps to enhance the adaptability of O. exigua to the environment (Bonin et al., 2007;Cen et al., 2023).The lack of significant differences among various geographical areas in O. exigua may result from the combination of life-history traits and ocean current circulation patterns (Leiva et al., 2023;Pérez-Portela et al., 2017;Taboada & Pérez-Portela, 2016).One key factor is the pelagic larval duration (PLD), which has a significant impact on gene flow and genetic structure in benthic organisms with a larval stage (Scheltema, 1986(Scheltema, , 1995;;Selkoe & Toonen, 2011).The planktotrophic larva of O. exigua remains in the water column for approximately 1 month (Kitazawa et al., 2015), exhibiting a strong dispersal ability that allows it to cross geographical barriers, thus significantly enhancing gene exchange.Additionally, ocean currents play a crucial role in larval dispersal (Scheltema, 1986).enabling gene exchange (Chen et al., 2020;Li, Yang, et al., 2020;Ramírez Ayala, ;Yang et al., 2019)  These transgressions correspond to high sea-level periods during  ka, an interglacial period), MIS6 (130-191 ka, possibly influenced by regional tectonic activity), and MIS5e (130-115 ka, the Eemian Interglacial) (Bintanja et al., 2002;Wang et al., 2020;Yi et al., 2012Yi et al., , 2013)).During glacial maxima, with a sea level drop of 130-150 m, the Bohai Sea and Yellow Sea Shelf were completely exposed (Liu et al., 2007;Xie et al., 1995).Therefore, we have sound reason to infer that different geographical populations of O. exigua might coalesce together due to the geographical range contraction during the glacial period and then expand population range when favorable conditions emerged during interglacial periods (Xi et al., 2023).The periodic coalescence of populations could enhance historical gene flow among populations (Xue et al., 2014).Expansion events in this period were welldocumented for other dominant species in the same regions, such as Nibea albiflora Richardson, 1846 (170-85 ka) (Han et al., 2008), Ophiura sarsii vadicola Djakonov, 1954 (111.8-60.1 ka) (Li, Dong, et al., 2020), and Atrina pectinata Linnaeus, 1767 (ca.160 ka) (Xue et al., 2014).The genetic population patterns of those species were similar to O. exigua, with high genetic diversity and low population divergence.
neutrality tests were used to infer the historical demography of O. exigua populations.A TCS haplotype network(Clement et al., 2002) was constructed using POPART v.1.7(Leigh & Bryant, 2015) to determine the genealogical relationships between haplotypes, visualizing the overall structure and spatial pattern of genetic diversity from the haplotype data of O. exigua.The mismatch distribution of populations from different locations was estimated by DnaSP to test the frequency distribution of pairwise differences under past population change.The demographic parameter tau (τ) was determined using DnaSP software to estimate the time since regional population expansion.The equation t = 2μk/τ(Rogers & Harpending, 1992) was used to determine the population size change of O. exigua, where t is the number of years since population expansion, μ is the per site per year substitution rate, and k refers to the gene fragment size.Finally, to assess changes in the size of the O. exigua population over time, we plotted the Bayesian Skyline Plots (BSP)

A
total of 80 COI, 72 NAD4, 75 ITS2, and 70 18S sequences were successfully amplified for 80 individuals.The average sizes of COI, NAD4, ITS2, and 18S gene fragments were 584, 645, 586, and 731 bp, respectively.All four gene sequences showed a high similarity (>98.97%) to the complete mitochondrial genome sequence of O. exigua in GenBank (ON729296).The average frequencies of the four nucleotides in the COI and NAD4 gene regions for all samples of O. exigua were 54.54% (A + T), 45.46% (G + C), 61.70% (A + T),and 38.30% (G + C), respectively.The A + T content was higher than the G + C content, which was consistent with the characteristics of the mitochondrial base composition(Galaska et al., 2019;Scouras et al., 2004;Sun et al., 2023).
) showed the connection among the six O. exigua populations.The ITS2 network had the simplest topology, with the highest number of shared haplotypes (H1) located at the center, while other haplotypes differed by only one or two substitutions from this main haplotype.Based on coalescent theory, H1 may represent the ancestral haplotype of O. exigua.In contrast, the COI and NAD4 networks revealed profound diverse bush-like genealogies, with black circles indicating missing haplotypes, and numerous haplotypes broadly and evenly distributed.The haplotypes from each sampling location were homogeneously mixed in the networks, failing to show distinct geographical patterns or population differentiation.
For instance, northward currents such as the Yellow Sea Warm Current and the Liaonan Coastal Current facilitate the northward dispersal of planktonic larvae from the DL population.Meanwhile, southward currents such as the Bohai Bay Coastal Current, Liaodong Bay Coastal Current, and Lubei Coastal Current promote the southward migration of planktonic larvae from QHD and RC populations, thereby

F
The mismatch distributions for all individuals in Northern China Sea.(a) COI and (b) NAD4.This figure was generated using DnaSP v6.12.03.F I G U R E 4 BSP Analysis of COI sequence datasets for all individuals in Northern China Sea.The X-axis represents time in millions of years before the present.The Y-axis represents the effective population size.The solid line indicates the median population size, and the 95% interval of HPD is shown in blue.This was generated using BEAUti v 2 and BEAST v 2.6.3.theresults of haplotype diversity, the neutral theory test(Tajima's D   and Fu's Fs), mismatch distribution and BSP all support the occurrence of a rapid population expansion event for O. exigua in the northern China Sea.The expansion time were estimated using different methods in this study.Mismatch analysis estimated the expansion time to be around 193,327 years ago, while BSP analysis dated the expansion event to 0.2-0.1 Ma.This timing closely aligns with three significant marine transgressions in the studied area since the Pleistocene.
, and other northwestern Pacific areas.Although O. exigua is widely distributed, this study focused solely on the northern China Sea.To unravel the population genetic patterns and evolutionary history of O. exigua in the northwestern Pacific areas, more samples from disparate geographical locations, especially the Japan Sea and the Korean Peninsula, need to be collected.The southern China Sea was not considered in our study because the southern specimens of this species in earlier records may have been misidentified.We confirmed significant morphological differences between our fresh O. exigua specimens and museum samples from southern China Sea (Marine Biological Museum, Chinese Academy of Sciences).Compared with our O. exigua specimens, the southern samples identified in Liao (2004) exhibit more sparsely distributed and elongated clubshaped spines on the disk, less densely covered radial shields with clearly visible outlines, more slender and transparent arm spines, and distinctive longitudinal black and white stripes on their dorsal arm plates.Based on morphological characteristics, the southern specimens should be assigned to Ophiothrix (Ophiothrix) ciliaris Genetic diversity statistics of Ophiothrix (Ophiothrix) exigua populations.

variation Df Sum of squares Variance components
Source of TA B L E 4 AMOVA analysis for the population of Ophiothrix (Ophiothrix) exigua.
. Frequent gene flow leads to the homogenization of genetic resources, resulting in the absence of significant genetic differentiation and typical geographical structure among O. exigua populations in the northern China Sea.
The population homogeneity found in O. exigua could be also the result of historical recolonization events, including habitat expansion and demographic growth during the Pleistocene period.In this study,